home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_10 / dirich.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  1.2 KB  |  45 lines  |  [MATF/MATL]
  1. function U = dirich(f1,f2,f3,f4,a,b,h,tol,max1)
  2. % U = dirich(f1,f2,f3,f4,a,b,h,tol,max1)
  3. % Dirichlet solution to Laplace's equation.
  4. % f1  is a boundary function, input.
  5. % f2  is a boundary function, input.
  6. % f3  is a boundary function, input.
  7. % f4  is a boundary unction, input.
  8. % a   is the width of [0 a], input.
  9. % b   is the width of [0 b], input.
  10. % h   is the step size, input.
  11. % tol is the tolerance, input.
  12. % max1 is the maximum number of iterations, input.
  13. % U  is the solution matrix, output.
  14. n = fix(a/h)+1;
  15. m = fix(b/h)+1;
  16. ave = (a*(feval(f1,0)+feval(f2,0)) ...
  17.     +  b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b);
  18. U = ave*ones(n,m);
  19. for j=1:m,
  20.   U(1,j) = feval(f3,h*(j-1));
  21.   U(n,j) = feval(f4,h*(j-1));
  22. end
  23. for i=1:n,
  24.   U(i,1) = feval(f1,h*(i-1));
  25.   U(i,m) = feval(f2,h*(i-1));
  26. end
  27. U(1,1) = (U(1,2) + U(2,1))/2;
  28. U(1,m) = (U(1,m-1) + U(2,m))/2;
  29. U(n,1) = (U(n-1,1) + U(n,2))/2;
  30. U(n,m) = (U(n-1,m) + U(n,m-1))/2;
  31. w = 4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2));
  32. err = 1;
  33. cnt = 0;
  34. while ((err>tol)&(cnt<=max1))
  35.   err = 0;
  36.   for j=2:(m-1),
  37.     for i=2:(n-1),
  38.       relx = w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+ U(i-1,j)-4*U(i,j))/4;
  39.       U(i,j) = U(i,j) + relx;
  40.       if (err<=abs(relx)), err=abs(relx); end
  41.     end
  42.   end
  43.   cnt = cnt+1;
  44. end
  45.